A Combinatorial Optimal Control Problem for Spacecraft Formation 
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Abstract — We consider a spacecraft formation reconfigura- 
tion problem in tlie case of identical spacecraft. This introduces 
in the optimal reconfiguration problem a permutational degree 
of freedom, in addition to the choice of individual spacecraft 
trajectories. We approach this using a coupled combinatorial 
and continuous optimization framework, in which the inner 
loop consists of computing the costs associated with a particular 
assignment by using a geometrically exact and numerically 
efficient discrete optimal control method based on Lie group 
variational integrators. In the outer optimization loop, combina- 
torial techniques are used to determine the optimal assignments 
based on the costs computed in the inner loop. The proposed 
method is demonstrated on the optimal reconfiguration problem 
for 5 identical spacecraft to go from an inline configuration to 
one equally spaced on a circle. 

I. Introduction 

The optimal control of spacecraft formations has received 
increased interest due to the NASA Terrestrial Planet Finder 
Project and the ESA Darwin Project. The objective is to use 
multiple spacecraft for cooperative missions such as long 
base-line interferometers, magnetosphere studies, and space- 
based communication networks. To accomplish various goals 
efficiently, it is often required to reconfigurate a formation 
during a mission. Since each spacecraft has limited fuel, it is 
important that the formation reconfiguration maneuvers are 
achieved with minimum fuel expenditure. 

Formation reconfiguration can be classified into two types: 
each spacecraft is required to be transferred into a specified 
location in the desired reconfigured formation, while in the 
other case, a specified location in the desired formation can 
be occupied by any spacecraft of a particular type [1]. In 
general, a formation is composed of identical spacecraft or 
groups of spacecraft at the same type, and the total fuel 
consumption depends on the permutations of the formation 
reconfigurations as well as the maneuver of each spacecraft 
to a specified location. 

In this paper, we study an optimal spacecraft formation 
problem integrated with an integer/combinatorial optimiza- 
tion approach for the assignment. Usually in combinatorial 
optimization problems for multiple agents, the dynamics of 
each agent is either ignored or simplified into an analytic 
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model such as a kinematics equation or a double integra- 
tor [2]. Here, we assume that each spacecraft evolves on the 
special Euclidean group SE(3), including both translational 
dynamics and rotational attitude dynamics under a central 
gravitational potential. Thus, finding the optimal control 
forces and moments on a spacecraft assigned to a fixed 
desired location is demanding even if the combinatorial 
assignment optimization problem is not considered. This 
is an interesting and challenging problem since it requires 
combining an integer/combinatorial optimization approach 
and an optimal control method over the non-trivial dynamics 
of spacecraft on SE(3). 

There has been some work in the literature on combina- 
torial optimization for spacecraft formation. The costs for 
all possible assignments are directly compared for space- 
craft moving along a straight path in [1]. This requires a 
large computational effort since optimal control problems 
associated with n! assignments have to be solved for a 
formation of n spacecraft. In [3], the special structure of 
a Hamiltonian system is utilized to expedite finding the 
solutions of optimization problems with varying boundary 
conditions. But, this requires a solution of the Hamilton- 
Jacobi partial differential equation. A stochastic optimization 
technique is used in [4]. 

To approximate the cost matrix used in the combinatorial 
assignment problem, we use the cost entries which have 
been explicitly computed, and the sensitivities of the cost 
to construct approximations to the remaining entries. The 
solution of the optimal control problem for each spacecraft is 
based on a structure-preserving numerical integrator referred 
to as a Lie group variational integrator [5]. Combined with 
an indirect optimization method, the Lie group variational 
integrator provides a geometrically exact but numerically 
efficient numerical optimization method for the dynamics 
of a rigid body [6]. The combinatorial optimization scheme 
for spacecraft formations that we present has the follow- 
ing important features: (1) dynamics of each spacecraft is 
nontrivial, (2) a discrete combinatorial optimization on a 
permutation group is explicitly integrated into the continuous 
optimal control problems, and (3) the problem is formulated 
and solved in a discrete time space using a Lie group 
variational integrator for overall computational accuracy and 
efficiency. 

This paper is organized as follows. Computational ap- 
proaches to solve an optimal control problem for a single 



spacecraft are summarized in Section and |III] Based 
on these results, a combinatorial optimization approach is 
developed in Section |IV] which is followed by a numerical 
example in Section IV] 

II. Lie Group Variational Integrator 

The configuration space for the translational and rotational 
motion of a rigid body is the special Euclidean group, 
SE(3) = R3(s)SO(3). We identify the cotangent bundle 
T*SE(3) with SE(3) x se(3)* by left translation, and we 
identify se(3)* with by an isomorphism between 
and se(3), and the isomorphism between se(3) and 5e(3)* 
induced by the standard inner product on M^. We denote the 
attitude, position, angular momentum, and linear momentum 
of the rigid body by {R, x, H, 7) G T*SE(3). 

The continuous equations of motion are given by 



m ' 



7 = / + u-^, 
R = RS{n), 

u + nxii = M- 



(1) 

(2) 
(3) 
(4) 



where il G M'^ is the angular velocity, and uf,u"''- £ are 
the control force in the inertial frame and the control moment 
in the body fixed frame, respectively. The constant mass of 
the rigid body is to G M, and J G M.^^^^ denotes the moment 
of inertia, i.e. 11 — Jil. The map S{-) : M.'^ ^ so(3) is an 
isomorphism between so (3) and M.^ defined by the condition 
S{x)y = X X y for all x,y G M"^. 

We assume that the potential is dependent on the position 
and the attitude; U{-) : SE(3) 1-^ R. The corresponding force 
and the moment due to the potential are given by 



J ' 

OX 

M — ri X + r2 X Ur^ + rs X UrZ, 



(5) 
(6) 

where rijUr^ G M"^ are the ith row vector of R and 
respectively. 

Since the dynamics of a rigid body has the structure of 
a Lagrangian or Hamiltonian system, they are symplectic, 
momentum and energy preserving. These geometric features 
determine the qualitative behavior of the rigid body dynam- 
ics, and they can serve as a basis for theoretical study of 
rigid body dynamics. 

In contrast, the most common numerical integration meth- 
ods, including the widely used (non-symplectic) explicit 
Runge-Kutta schemes, preserve neither the Lie group struc- 
ture nor these geometric properties. In addition, standard 
Runge-Kutta methods fail to capture the energy dissipation 
of a controlled system accurately [7]. Additionally, if we 
integrate (O by a typical Runge-Kutta scheme, the quantity 
R^R inevitably drifts from the identity matrix as the sim- 
ulation time increases. It is often proposed to parameterize 
rotations by Euler angles or unit quaternions. However, Euler 
angles are not global expressions of the attitude since they 
have associated singularities. Unit quaternions do not exhibit 



singularities, but are constrained to lie on the unit three- 
sphere S^, and general numerical integration methods do not 
preserve the unit length constraint. Therefore, quaternions 
have the same numerical drift problem. Renormalizing the 
quaternion vector at each step tends to break the conser- 
vation properties. Furthermore, unit quaternions, which are 
diffeomorphic to SU(2), double cover S0(3). So there are 
inevitable ambiguities in expressing the attitude. 

In [5], Lie group variational integrators are constructed 
by explicitly adapting Lie group methods [8] to the discrete 
variational principle [7]. They have the desirable property 
that they are symplectic and momentum preserving, and 
they exhibit good energy behavior for an exponentially long 
time period. They also preserve the Euclidian Lie group 
structure without the use of local charts, reprojection, or 
constraints. These geometrically exact numerical integration 
methods yield highly efficient and accurate computational 
algorithms for rigid body dynamics, and avoid singularities 
and ambiguities. 

Using the results presented in [5], a Lie group variational 
integrator on SE(3) for equations ([T]i-(|4|i is given by 



Xk H 7fc 

TO 
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{h + 4), (7) 
7fe + ^(/fe + 4)+^(/fe+i+4+i). (8) 



hS{Uk + - (Mfe + u^)) = FkJd - JdF't 



k 7 



Hi. 



fc+1 



F^nk + ^FnMk + ul^)+^l{Mk, 



(9) 
(10) 

^fe+i) J 
(11) 



where the subscript k denotes the /c-th step for a fixed 
integration step size G M. The matrix Jd G 



p3x3 



IS 



a nonstandard moment of inertia matrix defined by Jd = 
■^tr[ J] l3x3~J- The matrix Fk G S0(3) denotes the relative 
attitude between adjacent integration steps. 

For given {Rk,Xk, 11^,7^) and control input, (|9]l is solved 
to find Fk- Then {Rk+i,Xk+i) are obtained by ( fTOl i and 
(|7]l. Using ^ and (|6]l, {fk+i, Mk+i) are computed, and 
they are used to find (11^+1,7^+1) by (fTTT l and®. This 
yields a map (i?fe,a;fe, life, 7fc) 1-^ (i?fe+i, Xfe+i, Ilfe+i, 7^+1), 
and this process is repeated. The only implicit part is Q. 
The actual computation of Fk is done in the Lie algebra 
so (3) of dimension 3, and the rotation matrices are updated 
by multiplication. As this approach does not involve the 
component-wise integration of the kinematics equation (O, 
there is no excessive computational burden. It can be shown 
that this integrator has second-order accuracy. 

One of the distinct features of the Lie group variational 
integrator is that it preserves both the symplectic property 
and the Lie group structure of the rigid body dynamics. 
As such, it exhibits substantially improved computational 
accuracy and efficiency compared with other geometric in- 
tegrators that preserve only one of these properties, that is 
symplectic Runge-Kutta methods that do not preserve Lie 
group structure or non-symplectic Lie group methods [9]. 



The symplectic property of numerical integrators is important 
even in the case of controlled dynamics, since the dissipation 
rate of the total energy is typically computed inaccurately by 
non-symplectic integrators [7]. 

III. Optimal Control of a Rigid Body on SE(3) 

We first summarize a computational approach to solve 
the optimal control problem for a single rigid body in 
which the translational dynamics and the rotational attitude 
dynamics are coupled [6]. This approach is extended to solve 
a combinatorial optimal formation problem for multiple rigid 
bodies in Section HV] 

A. Problem Formulation 

An optimal control problem is formulated for the 
maneuver of a rigid body from a given initial con- 
figuration {Rq,xo,IIq,jo) to a desired configuration 
{R%, x%, 11%, jff) during a given maneuver time N. Control 
inputs are parameterized by their value at each time step. The 
performance index is the square of the weighted I2 norm of 
the control inputs. 

given: (xq, 70, i?o, IIo), (a;^, 7d, ^d, n^), N, 

N-l . , , N 

such that {xn,1n,Rn,'^n) = {xd,ld,Rd,'n.d), 
subject to discrete equations of motion (l7b-(fTT]). 

where Wf , Wm G M"^ ^ are symmetric positive definite 
matrices. 

B. Computational Approach 

We solve this optimal control problem using an indirect 
method; necessary conditions for optimality are obtained by 
using variational expressions that respect the geometry of 
the configuration space, and the corresponding two point 
boundary value problem is solved by using the shooting 
method. Here we use a modified version of the discrete 
equations of motion with first order accuracy, because it 
yields a compact form for the necessary conditions. 

Define an augmented performance index as 

fc=0 

h 



rm„,m 



A 



l.T 



-Xfc+i + Xfc H 7fc 

m 



+ A 



2,T 



|-7fc+i +lk + hfk+1 + 



+ \l'^S-'{\ogmiFk-RlRk+i)) 

+ 4'^ {-Hfc+i + FjHfe + h (A/fe+i + } , 

where X^, G M'^ are Lagrange multipliers. The constraint Q 
is considered implicitly using a constrained variation. Setting 
SJ'a = for all variations, we obtain necessary conditions 
for optimality as follows. 



fc+iAfc+1, 



(13) 
(14) 



Afc = [Ai;A2;A|;A|.] e M^^^ and e Ri2xi2 

is suitably defined in terms of (i?^, x^, 11^, 7^). 
Together with the discrete equations of motion, 
this yields a map {{Rk, xo,Uk,jk), Xk} ^ 
{(i?fc+i, a::fc+i,nfc+i,7fc+i), Afe+i}. 

The necessary conditions for optimality are expressed in 
terms of a two point boundary problem on T*SE(3) and 
its dual. This problem is to find the optimal discrete flow, 
multiplier, and control inputs to satisfy the equations of 
motion, optimality conditions (fT2]i.(fT3Tl, multiplier equations 
(fT4l l. and boundary conditions simultaneously. 

We use the shooting method [10]. A nominal solution 
satisfying all of the necessary conditions except the boundary 
conditions is chosen. The unspecified initial multiplier is up- 
dated by successive linearization so as to satisfy the specified 
terminal boundary conditions in the limit. The optimality 
conditions (fT2] l and (flJl i are substituted into the equations 
of motion and the multiplier equations. The sensitivities of 
the specified terminal boundary conditions with respect to 
the unspecified initial multiplier conditions is obtained by a 
linear analysis. 

Let Zk e be the variation of the state given by Zk = 
[£,k',^Xk;SIlk]S^k], where (k G denotes the variation 
of the rotation matrix as SRk = ^_^Rk exp S (Ck) = 
RkS{Ck), and variations for other variables are defined in 
the usual sense. The linearized equations of motion and the 
linearized multiplier equation can be written as 



A'^SXk, 



Zk+i = AkZk ■ 



5Xk — Ak+iZk+i + Aj^^i5Xk+i 



(15) 
(16) 



where Al^,Al\i G M^^^^^ can be computed expHcitly. The 
solution of the linear equations ( fTSl l and ( fTSI l can be obtained 



Zk 

SXk 



SXa 



(17) 



where <^>'^ G R^^^^^ 

For the given two point boundary value problem, zq = 
since the initial condition is fixed, and Aat is free. Thus, 



ZN = '^nSXq. 



(18) 



W^'Xi, 



(12) 



The matrix represents the sensitivity of the specified 
terminal boundary conditions with respect to the unspecified 
initial multipliers. Using this sensitivity, an initial guess 
of the unspecified initial conditions is iterated to satisfy 
the specified terminal conditions in the limit. Any type 
of Newton iteration can be applied. We use a line search 
with backtracking algorithm, referred to as Newton-Armijo 
iteration in [11]: the outer loop computes the sensitivity 
derivatives to obtain the Newton search direction, and the 
inner loop performs a line search to find the largest step size 
along the given search direction. 



C. Properties of Computational Approach 

The key feature of this computational approach for the 
optimal control problem of a single rigid body is that it is 
discretized from the problem definition level using the Lie 
group variational integrator This is in contrast to obtaining 
continuous time necessary conditions, which are discretized 
to numerically solve the two point boundary value prob- 
lem. In this computational approach for the optimal control 
problem, the discrete necessary conditions for optimality are 
obtained by a variational principle. 

The main advantage of the shooting method is that the 
number of iteration variables, the initial Lagrange multi- 
pliers, is small. In other approaches, an initial guess of a 
control input history or multiplier history are iterated, so 
the number of optimization parameters is proportional to 
the number of discrete time steps. The difficulty is that 
the extremal solutions are sensitive to small changes in the 
unspecified initial multiplier values. The nonlinearities also 
make it hard to construct an accurate estimate of sensitivity, 
perhaps resulting in numerical ill-conditioning. 

Here, the discrete necessary conditions for optimality 
preserve the geometric structure of the optimal control 
problem. Thus, there is no geometrical error introduced 
by the numerical integration algorithm itself. It turns out 
that, combined with the shooting method, this computational 
approach provides a geometrically exact and numerically 
efficient solution to this highly nonlinear, non-convex rigid 
body optimal control problem [6], [12]. This is used as a 
basic tool for the combinatorial optimization problem for 
multiple rigid bodies. 

IV. Optimal Formation Control of Rigid Bodies 
A. Problem Formulation 

We study an optimal formation control problem of n 
identical rigid bodies where the maneuver of each body 
is described by (l7])-(fTT]). The objective is to find the opti- 
mal control forces and moments for each rigid body such 
that the group moves from a given initial configuration 
(i?Q,a;g,nQ,7g) for i G {l,2,...,n} to a desired target 
T G K''" during a given maneuver time N, where the 
superscript i denotes the i-th rigid body. 

More precisely, we assume that the n desired positions 
which all rigid bodies are located at the 
terminal maneuver time, are given as functions of parameters 
G M'. The desired attitude, the linear momentum, and the 
angular momentum at the terminal time, (_Rrf,nd,7d), are 
assumed to be fixed and to be the same for all rigid bodies. 

Since all rigid bodies are identical, there are n\ possible 
combinatorial assignments for n rigid bodies to these n 
desired locations. Let {a^} be a n x n matrix composed 
of binary elements {0, 1}, referred to as an assignment or a 
permutation matrix. Each element of the assignment matrix 
aij represents the assignment of the i-th rigid body to the j- 
th desired terminal position a;;^. If = 1, the i-th rigid body 
is assigned to the j-th node, and if — 0, the i-th rigid 
body is not assigned to the j-th node. Thus, the assignment 



is valid if X]J=i '^y — SILi o^ij = 1- The assignment matrix 
can be equivalently expressed as a set A = {(z, j)|ay = 1}, 
and the particular desired points assigned by the i-th rigid 
body for an assignment A is denoted by Ai G {1, . . . ,n}. 
In other words, for an assignment A, the i-th rigid body is 
assigned to the Ai-i\\ desired location, x^' . 

The target is defined in terms of a parameter 6 and an 
assignment A as follows. 

T{e,A) = [xf{e)y_^ gm^". 

Thus, for a given parameter G M' and a given assignment 
A G Sn, the terminal boundary conditions for all rigid bodies 
are completely determined. 

The performance index is the sum of the squares of the 
weighted I2 norms of the control inputs. The optimal control 
problem for a formation of n rigid bodies is formulated as 

given: {(x^, 7^, i?j„ nj,)};^^ , , i?,, H^, 7<i), iV, 

i=l k=0 

such that ^^{R%,x%n%,j%) = {Rd,x^'{0),Udnd)] , 
subject to discrete equations of motion (l7b-(fTT]). 

where Wf,Wm S M'^^'^ are symmetric positive definite 
matrices. 

Since we have neglected the gravitational interactions 
between the rigid bodies, the dynamics of the rigid bodies 
are only coupled through the terminal boundary conditions. If 
the parameter 9 and the assignment A are prescribed so that 
the terminal configurations of all rigid bodies are completely 
determined, then the optimal control problems for n rigid 
bodies can be solved independently using the computational 
approach presented in Section |III] The formation cost is the 
summation of the resulting costs of each rigid body. There- 
fore, the optimal formation control problem for multiple rigid 
bodies consists in finding the optimal value of the parameter 
and the optimal assignment of rigid bodies among the n\ 
possible assignments. This problem formulation is similar 
to the optimal formation reconfiguration problem presented 
in [13] except that we include the combinatorial assignment 
problem explicitly in this paper 

This requires combining the optimal control approach and 
the integer/combinatorial assignment optimization over the 
non-trivial dynamics of rigid bodies on SE(3)". We present 
a computational approach to this integrated optimal control 
problem. 

B. Optimal Control of Rigid Bodies on SE(3)" 

We first solve the optimal formation control problem 
assuming that an assignment A is pre-determined and fixed. 
We use a hierarchical optimal control approach [13]. Since 
the parameter 9 completely defines the terminal configuration 
of all rigid bodies for the fixed assignment A, it also 
determines the corresponding cost by taking the sum of 
the cost of the optimal trajectories for each rigid body. 



Thus, the optimization problem is decomposed into an outer 
optimization problem to find the optimal value of 9 that 
minimizes the total cost, and an inner optimization problem 
to find the optimal control forces and moments for the given 
value of 6. This is a consequence of the fact that, 

min J(w,6l) = min(min{J(M,6l)|6l = e'}\ . (19) 
u,e 9' K u J 

The inner optimization problem is solved by using the 
computational approach given in Section The optimal 
value of 6 is found by using a parameter optimization method 
with an explicit expression for the gradient. 

Sensitivity of the cost: Based on the solution of the 
optimal control problem of the i-th body, the sensitivity of 
the cost with respect to the parameter can be obtained as 
follows. Let £ R be the contribution of the i-th body to 
the performance index so that J = X]"=i 



E 

fc=0 



\iuf:l^Wf4:l, + ^^^w^u^-,. (20) 



Suppose that the variation of the initial condition and the ter- 
minal boundary condition are given by zl^^z)^ S M^^. Using 
( [TtI i. the corresponding variation of the initial multiplier for 
an optimized solution is given by 



5\\ 



'^,^{<^'^^T\-^n''4 + ^n)- (21) 

Substituting this into ( [TtT i. we obtain the variation of the 
multiplier as 



21,i. 



). (22) 



Since the control input is expressed in terms of the multiplier 
by the optimality condition (fT2]i.(fT3]), the costs given in (l20b 
are represented as follows. 



N-l 
fc=0 



4,i\ 



wx 



fc=0 



where 



- '^-^k '--^ 



diag[03x3,M^7\03x3,W'-^] e 
variation of the cost is given by 



Tn>12xl2 



and W = 
Using dill, the 



N-l 

k=0 
N-l 



22,1 



+ 



k=0 
N-l 

h J2 iKfw<i> 

I k=0 



N 



(23) 



These represent the sensitivities of the optimal costs with 
respect to the initial condition and the terminal boundary 
condition; gK^^^i^ 

' azi, ' dz' 



Computational approach: The sensitivity of the perfor- 
mance index with respect to the target space parameter is 
given by 



dj__^dc_ 



89 



09 



(24) 



where e M^^'^ is composed of the fourth to sixth 



dx 

elements of 

dzl, 



1 dc^ 

IT 10)12x3 



4^B, where B 

dz' 



[03x3,-^3x3,03x3,03x3]"' G 

A quasi-Newton method can be applied to the outer 
optimization using this gradient computation. For example, 
the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method to 
solve an unconstrained nonlinear optimization problem with 
an approximated Hessian is summarized as follows [11]. 



1 

2 
3 
4 
5 
6 
7 
8 
9 
10 



Guess an initial parameter 9. 
Solve the n optimal control problems. 
Find 4# using ( |24] |. 



H 



while ll^ll > e 

Find a line search direction; D = 
Perform line search 6 = 6 + aD 
Update the Hessian H. 

Solve the n optimal control problems. 

aj 

06 



-Id J 



Find ^ using ( |24] |. 



end while 



Here e denotes a stopping criterion and a is a scaling factor, 
respectively. The Hessian can be initialized with H = J^xi 
and updated during the BFGS iterations. 

The major computational burden is in the third step and the 
fifth step. The computation time for the numerical iteration 
presented in Section |lll] can be substantially reduced if 
we have a good guess of the initial multiplier. At each 
iteration, we store the optimized initial multiplier for the 
terminal boundary conditions, and we use the accumulated 
data to initialize the initial multiplier at the next iteration. For 
example, UaS can be used to find an educated guess of the 
initial multiplier for neighboring boundary conditions. This 
reduces the computational burden as the iterations proceed, 
which will be shown by a numerical example in Section [V] 

C. Assignment Optimization Problem 

Now, we solve the optimal formation control problem 
assuming that the target parameter 9 is determined and fixed. 
For the given value of 9, the n desired points \^^d\ i-v 
which all rigid bodies are located at the terminal maneuver 
time, are completely defined. Thus, there are n\ possible 
combinatorial assignments. 

Let {c*-' } be a n X n matrix, referred to as a cost matrix. 
Each element c*^ represents the optimal cost of the i-th 
rigid body transferred to the j-th desired location. For an 
assignment {a^j}, the performance index is given by J7 = 
^YTi ^=1 c*"'aii- The optimal assignment problem is given by 



Subject to ^ Uij = 1 for i 6 {1, 2, . . . , n} . 



^ay = l for j G {1,2, . . . , 

1=1 

£{0,1} fori,j e {l,2,...,n}. 

Since we assume that there is no interaction between rigid 
bodies, the cost matrix is independent of the assignment. 
For the given value of the target parameter 9, we must 
solve at most optimal control problems to obtain the 
cost matrix. Once we have the complete cost matrix, the 
optimal assignment can be obtained by comparing costs for 
all possible assignments or by using the Hungarian method 
for large dimensional systems [14]. 

It is often expensive to obtain the cost matrix. Each 
element of the cost matrix is a solution of the optimal 
control problem presented in Section For the formation 
optimization problem, we need to find the cost matrix with 
varying values of the target parameter 9. Thus, the objective 
of this subsection is to find the optimal assignment without 
solving all optimal control problems. We start with an 
initial single spacecraft optimal trajectory computation, and 
use its optimal cost and the sensitivities given in ( l23T l to 
populate the remaining entries of the cost matrix. 

Suppose that we solve the optimal control problem of the 
first rigid body transferred to the first desired point to obtain 
c^^. Since this optimal solution is obtained by computing 
the linearized equations in iT7[ . we can find the sensitivity 
of the cost with respect to the terminal boundary condition 
by using (|23] i. without need of additional computational 
burden. Then, the optimal cost of transferring the first body 
to the other desired points, say Xd, is approximated as 



TABLE I 

Methods to choose sensitivities 



dc 



1 



2„11 



;Axd, (25) 



dxd 2 

l^. The first order sensitivity 



da 



T- IS 



where Ax^; = Xd — x^ 

is computed exactly, and the second order Hessian 
initially set to zero, and the approximation is improved as 
other optimal solutions become available. For example, if the 
solution of the optimal control problem of the first rigid body 
transferred to the second desired point is found, we obtain 
the exact value of c^^ and . This provides the following 
1 and 3 dimensional constraints on the Hessian, 



dxd 



Thus, the 6 elements of the Hessian can be approximated in 
either the minimum norm or least squares sense if additional 
optimal solutions involving the first rigid body are available. 

This approach approximates the elements of the cost ma- 
trix along rows using the sensitivity of the cost with respect 
to terminal boundary conditions. A similar approximation 
along columns can be made by using the sensitivity of the 
cost with respect to the initial conditions. In the combinato- 
rial optimization process, we utilize both approximations in 
order to avoid local minima. The advantage is that we use 



Method 


Sensitivity Selection 


Term. 


Use terminal sensitivity always. 


Init. 


Use initial sensitivity always. 


Rand. 


Choose one of sensitivities randomly. 


Rpt 


If the same assignment is repeated two times, switch to the 
different sensitivity. 


Alt. 


Alternate sensitivities. 


Comp. 


For each element of the cost matrix, compare the number of 
available solutions along the row direction, and the number 
of available solutions along the column direction. Select one 
of sensitivities that has more available solutions. If both di- 
rections have the same number of available solutions, choose 
randomly. 



all of the sensitivity information available up to the current 
iteration in order to estimate the new cost matrix. 

We construct a combinatorial assignment method using 
these approximations. 

i) Guess an initial assignment and solve the correspond- 
ing optimal control problems for this assignment. 

ii) Estimate the cost matrix using the linear approxima- 
tions. 

iii) Find a new assignment using the estimated cost matrix, 
and solve the corresponding optimal control problems 
for this assignment. 

iv) Find the best assignment using all of the solutions of 
the optimal control problems obtained so far 

v) Construct a second-order approximation equation dZSl) 
based on the best assignment, and estimate the cost 
matrix. 

vi) Repeat (iii)-(v) until the same approximation is re- 
peated M times in a row at (v). 

Numerical simulations show that setting M = 3 is sufficient 
to find the global optimal assignment. At Steps (ii) and (v), 
we construct the approximation to the cost matrix using 
either the sensitivity of the cost with respect to the initial 
conditions or the sensitivity of the cost with respect to the 
terminal conditions. Several ways are summarized in Table 
U Numerical simulations demonstrate that using both types 
of sensitivities has advantages and the last method is more 
efficient than the others in terms of finding the global optimal 
assignment with minimal computational effort. 

D. Computational Approach for Optimal Formation Control 

We have presented two optimization approaches; finding 
the optimal value of the space parameter for a given assign- 
ment, and finding the optimal assignment for a given value 
of the target parameter. We integrate both methods using a 
hierarchical optimization approach similar to ( fT9] l. 

The original optimization problem is stated as finding 
the optimal control forces, moments, target parameter, and 
assignment that minimizes the total cost. 

min J(u,9,A). 

u,e,A 



Equivalently, this can be stated as finding the optimal assign- 
ment over the solutions for the optimal control inputs and 
the target parameters 

min \ min {J(u, e,A)\A^A'}\. 
A' u,e J 

In the inner stage, we optimize the target parameter and the 
control inputs using the continuous optimization approaches 
presented in Section IIV-BI and in the outer stage, we find 
the optimal assignment using the combinatorial optimiza- 
tion approach presented in Section IIV-CI The optimization 
process is terminated when the iterations yield a solution 
that is optimal for both the continuous and combinatorial 
optimization stages. 

V. Numerical Example 

We study a maneuver involving 5 identical rigid spacecraft 
under a central gravity field. We assume that the mass of each 
spacecraft is negligible compared to the mass of a central 
body, and we consider a fixed frame attached to the central 
body as an inertial frame. The resulting model is a Restricted 
Full Body Problem (RFBP) [6]. 

Each spacecraft is modeled as a dumbbell, which consists 
of two equal spheres and a massless rod. The gravitational 
potential is given by 

U .±-^±^-1—^^. (26, 

l—l q—l II ' y II 

where G G K is the gravitational constant, M, to G R are 
the mass of the central body, and the mass of the dumbbell, 
respectively. The vector G M.^ is the position of the qth 
sphere from the mass center of the i-th dumbbell expressed 
in the body fixed frame (q G {1, 2}). The mass, length, and 
time dimensions are normalized by the mass of the dumbbell, 
the radius of a reference circular orbit, and its orbital period. 

The spacecraft are initially aligned along a radial direction 
as shown in Fig. |l(a)| At the terminal time, we require that 
spacecraft are equally distributed on a target circle described 
by the location of its center Xo G M.^, the radius Tq G M, 
and the unit normal vector Uo G Let 6*' G be the 
angle of the i-th spacecraft on the target circle from a given 
reference direction as shown in Fig. |l(b)| We choose the 
target parameter as the angle of the first rigid body. The 
target T is given by 

T{9^,A) = {xo+^o cos6''ei-^ro sin6''e2}"^^ , 

where ei = -jr^, £2 — eiX Uo are unit vectors in the target 

ll^oll 

plane, and the angle 6'* is chosen to distribute the spacecraft 
uniformly on the circle 

e\9\A) = e^ + ^{A,-i). 
5 

Since the target parameter 9^ determines the terminal posi- 
tion of the first spacecraft completely, we require that the 
first spacecraft be assigned to the first desired location, i.e. 
Ai = 1. There remains 4! assignments for the other four 
spacecraft. Thus, the optimization parameters are the angle 




(a) The initial formation and the terminal formation 
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(b) The terminal formation on a target circle 



Fig. 1. The initial formation and the desired terminal formation of 5 
dumbbell spacecraft on a circle 



of the first spacecraft on the target circle, the 4! assignments 
for the remaining spacecraft, and the control inputs and 
moments. 

The iteration procedure for a particular numerical imple- 
mentation of the optimization are shown as follows. The 
target parameter, assignment, cost, and computation time on 
an Intel Pentium M 1.73GHz processor using MATLAB are 
given for each iteration step. 

i) The initial guess of the assignment is given by ^ = 
{(1,1), (2,4), (3,2), (4,3), (5,5)}. 
ii.a) For the given assignment, the optimal value of 9^ = 
2.4520 is obtained in 48.32 minutes with cost J = 
8.6984. 

ii.b) For the given value of 9^, the optimal assignment of 
A = {(1,1), (2, 5), (3, 2), (4, 3), (5, 4)} is obtained in 
3.04 minutes with cost J = 8.6905. 

ii. c) For the given 9^ and the given assignment, we check 

1^ = 2.42 X 10"^. Repeat iteration. 

iii. a) For the given assignment, the optimal value of 9^ = 

2.5084 is obtained in 12.69 minutes with cost J' = 
8.6898. 

iii.b) For the given value of 9^, the same optimal assignment 
of A = {(1, 1), (2, 5), (3, 2), (4, 3), (5, 4)} is obtained 
in 2.98 minutes with cost J = 8.6898. 

iii.c) For the given 9^ and the corresponding optimal assign- 
ment, we check |gT = 9.99 x 10^^. 
iv) The optimization is terminated with optimal cost 
J = 8.6898 for 9^ = 2.5084 and A = 
{(1, 1), (2, 5), (3, 2), (4, 3), (5, 4)} in total computation 




Fig. 2. Optimal spacecraft formation reconfiguration maneuver 

time 67.03 minutes. 
The corresponding maneuvers for all the spacecraft are 
shown in Fig. |2] The computation time for optimizing the 
target space parameter is reduced from 48.32 minutes at Step 
(ii.a) to 12.69 minutes at Step (iii.a). At each iteration, we use 
the optimization data accumulated in the previous iterations 
in order to initialize the initial multiplier for the optimal 
control problems. This reduces the computation time as the 
iterations proceed. 

In order to estimate the distribution of the possible so- 
lutions, we uniformly discretize the interval [0, 27r) by 100 
points for the target parameters and we find the total costs of 
4! assignments for each value of the target parameter The 
histogram for total costs of the corresponding 100 x 4! = 
2400 solutions is shown in Fig. |3(a)| 

Numerical simulations show that the numerical optimized 
solution obtained depends on the initial guess of the as- 
signment, and it is independent of the initial guesses for 
the target parameter and the initial Lagrange multiplier We 
repeat the numerical optimization for all possible 4! initial 
guesses of the assignment. Fig. |3(b)| shows the histogram of 
the optimized total costs for varying initial guesses for the 
assignment. Six initial assignments converged to the global 
optimal solution with J = 8.6898. 

VI. Conclusions 

A combinatorial optimization method for spacecraft for- 
mation flight is presented. The objective is to transfer a group 
of identical spacecraft to a desired formation with mini- 
mum fuel expenditure. The assignment optimization over the 
discrete permutation group is explicitly integrated with the 
solutions of optimal control problems for combined orbital 
and rotational maneuvers of spacecraft, which are described 
by a Lie group variational integrator The computational 
efficiency of the presented method is demonstrated by means 
of a numerical example. 
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